// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

#include <avtPolyhedralSplit.h>
#include <vtkDataArray.h>
#include <vtkUnsignedIntArray.h>
#include <map>

#include <DebugStream.h>

// ****************************************************************************
// Method: avtPolyhedralSplit::avtPolyhedralSplit
//
// Purpose:
//   Constructor
//
// Arguments:
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:16:41 PST 2010
//
// Modifications:
//   Brad Whitlock, Tue Oct 26 16:27:26 PDT 2010
//   I removed the arguments and made polyhedralSplit into a vector.
//
// ****************************************************************************

avtPolyhedralSplit::avtPolyhedralSplit() : polyhedralSplit(), nodesForPolyhedralCells()
{
}

// ****************************************************************************
// Method: avtPolyhedralSplit::~avtPolyhedralSplit
//
// Purpose:
//   Destructor.
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:17:13 PST 2010
//
// Modifications:
//
// ****************************************************************************

avtPolyhedralSplit::~avtPolyhedralSplit()
{
}

// ****************************************************************************
// Method: avtPolyhedralSplit::Destruct
//
// Purpose:
//   Static destructor function
//
// Arguments:
//   ptr : Pointer to the avtPolyhedralSplit object to delete.
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:17:27 PST 2010
//
// Modifications:
//
// ****************************************************************************

void
avtPolyhedralSplit::Destruct(void *ptr)
{
    avtPolyhedralSplit *This = (avtPolyhedralSplit *)ptr;
    delete This;
}

// ****************************************************************************
// Method: avtPolyhedralSplit::AppendCellSplits
//
// Purpose:
//   Record the cellid and number of splits for the next polyhedral cell.
//
// Arguments:
//   cellid  : The original cell id of the polyhedral cell.
//   nsplits : The number of splits for the polyhedral cell.
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:17:55 PST 2010
//
// Modifications:
//
// ****************************************************************************

void
avtPolyhedralSplit::AppendCellSplits(int cellid, int nsplits)
{
    polyhedralSplit.push_back(cellid);
    polyhedralSplit.push_back(nsplits);
}

// ****************************************************************************
// Method: avtPolyhedralSplit::AppendPolyhedralNode
//
// Purpose:
//   Store a polyhedral node id.
//
// Arguments:
//   id : The id of the node to save off.
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 16:50:25 PST 2010
//
// Modifications:
//
// ****************************************************************************

void
avtPolyhedralSplit::AppendPolyhedralNode(int id)
{
    nodesForPolyhedralCells.push_back(id);
}

// ****************************************************************************
// Method: avtPolyhedralSplit::ExpandDataArray
//
// Purpose:
//   Return a copy of the input data array with polyhedral cells' data
//   replicated the appropriate number of times.
//
// Arguments:
//   input        : The input data array.
//   zoneCent     : True if the data is zone centered.
//   averageNodes : If true then average the node data for new polyhedral nodes.
//                  If false, use the dominant value.
//
// Returns:    A new data array
//
// Note:       This just works for zonal data right now.
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:19:07 PST 2010
//
// Modifications:
//   Brad Whitlock, Tue Oct 26 16:23:45 PDT 2010
//   I turned polyhedralSplit into a vector. Allow for the input to be returned
//   in the case where there are now polyhedral cells. I also added a new case
//   for figuring out the new nodal values for new polyhedral nodes so we can
//   use the dominant value. This lets us handle nodal material data.
//
//   Brad Whitlock, Mon Nov 22 14:14:16 PST 2010
//   Set the component count so it works on vector data too.
//
// ****************************************************************************

vtkDataArray *
avtPolyhedralSplit::ExpandDataArray(vtkDataArray *input, bool zoneCent,
    bool averageNodes) const
{
    const char *mName = "avtPolyhedralSplit::ExpandDataArray: ";
    vtkDataArray *output = NULL;
    int polyhedralCellCount = (int)polyhedralSplit.size() / 2;

    if(polyhedralCellCount == 0)
    {
        debug4 << mName << "Returning input data array " << endl;
        output = input;
        output->Register(NULL);
    }
    else if(zoneCent)
    {
        debug4 << mName << "Expanding zonal data" << endl;
        output = input->NewInstance();
        int bloat = 0;

        for(size_t i = 1; i < polyhedralSplit.size(); i += 2)
           bloat += (polyhedralSplit[i] - 1);
        output->SetNumberOfComponents(input->GetNumberOfComponents());
        output->SetNumberOfTuples(input->GetNumberOfTuples() + bloat);
        output->SetName(input->GetName());

        vtkIdType out = 0;
        int phIndex = 0;
        for(vtkIdType cellid = 0; cellid < input->GetNumberOfTuples(); ++cellid)
        {
            if(phIndex >= polyhedralCellCount ||
               cellid < polyhedralSplit[phIndex*2])
            {
                output->SetTuple(out++, input->GetTuple(cellid));
            }
            else
            {
                int nrepeats = polyhedralSplit[phIndex*2+1];
                for(int j = 0; j < nrepeats; ++j)
                    output->SetTuple(out++, input->GetTuple(cellid));
                phIndex++;
            }
        }
    }
    else
    {
        output = input->NewInstance();
        output->SetNumberOfComponents(input->GetNumberOfComponents());
        output->SetNumberOfTuples(input->GetNumberOfTuples() + polyhedralCellCount);
        output->SetName(input->GetName());
        // Copy all of the original node data
        vtkIdType out = 0;
        for(vtkIdType nodeid = 0; nodeid < input->GetNumberOfTuples(); ++nodeid)
            output->SetTuple(out++, input->GetTuple(nodeid));

        // Create average values for the new nodes.
        if(averageNodes ||
           input->GetNumberOfComponents()>1/*For ease*/)
        {
            debug4 << mName << "Expanding nodal data using averaging for the "
                "new node values" << endl;
            double *avg = new double[input->GetNumberOfComponents()];
            for(size_t i = 0; i < nodesForPolyhedralCells.size(); )
            {
                int nnodes = nodesForPolyhedralCells[i++];
                memset(avg, 0, input->GetNumberOfComponents() * sizeof(double));
                for(int j = 0; j < nnodes; ++j)
                {
                    int nodeid = nodesForPolyhedralCells[i++];
                    double *val = input->GetTuple(nodeid);
                    for(int c = 0; c < input->GetNumberOfComponents(); ++c)
                        avg[c] += val[c];
                }
                for(int c = 0; c < input->GetNumberOfComponents(); ++c)
                    avg[c] /= double(nnodes);
                output->SetTuple(out++, avg);
            }
            delete [] avg;
        }
        else
        {
            debug4 << mName << "Expanding nodal data using dominant value for "
                "the new node values" << endl;

            // Just 1-component values for now.
            for(size_t i = 0; i < nodesForPolyhedralCells.size(); )
            {
                int nnodes = nodesForPolyhedralCells[i++];
                // Count up the unique values.
                std::map<double,int> counts;
                std::map<double,int>::iterator it;
                for(int j = 0; j < nnodes; ++j)
                {
                    int nodeid = nodesForPolyhedralCells[i++];
                    double val = input->GetTuple1(nodeid);
                    it = counts.find(val);
                    if(it == counts.end())
                        counts[val] = 1;
                    else
                        counts[val]++;
                }

                // Figure out which is the most common value.
                it = counts.begin();
                int maxcount = it->second;
                double mostCommon = it->first;
                for(; it != counts.end(); ++it)
                {
                    if(it->second > maxcount)
                    {
                        maxcount = it->second;
                        mostCommon = it->first;
                    }
                }

                // Save the most common value as the new node value
                output->SetTuple1(out++, mostCommon);
            }
        }
    }

    debug4 << mName << "Output data has " << output->GetNumberOfTuples()
           << " tuples" << endl;

    return output;
}

// ****************************************************************************
// Method: avtPolyhedralSplit::CreateOriginalCells
//
// Purpose:
//   Creates an original cells array based on the polyhedral split data.
//
// Arguments:
//   domain          : The current domain number.
//   normalCellCount : The number of normal cells.
//
// Returns:    An original cells array.
//
// Note:
//
// Programmer: Brad Whitlock
// Creation:   Fri Mar 12 15:20:24 PST 2010
//
// Modifications:
//   Brad Whitlock, Tue Oct 26 16:23:45 PDT 2010
//   I turned polyhedralSplit into a vector.
//
// ****************************************************************************

vtkDataArray *
avtPolyhedralSplit::CreateOriginalCells(int domain, int normalCellCount) const
{
    vtkUnsignedIntArray *originalCells = vtkUnsignedIntArray::New();
    originalCells->SetNumberOfComponents(2);
    int bloat = 0;
    int polyhedralCellCount = (int)polyhedralSplit.size()/2;
    for(int i = 0; i < polyhedralCellCount; ++i)
        bloat += polyhedralSplit[i*2+1];
    originalCells->SetNumberOfTuples(normalCellCount + bloat);
    originalCells->SetName("avtOriginalCellNumbers");

    int phIndex = 0;
    unsigned int *oc = (unsigned int *)originalCells->GetVoidPointer(0);
    int allCells = normalCellCount + polyhedralCellCount;
    for(int origCell = 0; origCell < allCells; ++origCell)
    {
        oc[1] = origCell;
        if(phIndex >= polyhedralCellCount)
        {
            *oc++ = domain;
            *oc++ = origCell;
        }
        else if(origCell < polyhedralSplit[phIndex*2])
        {
            *oc++ = domain;
            *oc++ = origCell;
        }
        else
        {
            int nrepeats = polyhedralSplit[phIndex*2+1];
            for(int j = 0; j < nrepeats; ++j)
            {
                *oc++ = domain;
                *oc++ = origCell;
            }
            phIndex++;
        }
    }

    return originalCells;
}

